The time-dependent Gutzwiller theory for 

multi-band Hubbard models 

b' 

(N 
sj^ , E. V. Oelsen, G. Seibold, and J. Biinemann 

O ' Institut ffir Physik, BTU Cottbus, P.O. Box 101344, 03013 Cottbus, Germany 

OQ ' Abstract. We formulate a multi-band generalisation of the time-dependent 

^^ . Gutzwiller theory. This approach allows for the calculation of general two-particle 

response functions, which are crucial for an understanding of various experiments in 
Q ' solid-state physics. As a first application, we study the momentum- and frequency- 

j" , resolved magnetic susceptibility in a two-band Hubbard model. Like in the underlying 

"ti ' ground-state approaches we find significant differences between the results of our 

method and those from a time-dependent Hartree-Fock approximation. 

Td ■ PACS numbers: 71.10.-w,71.10.Fd,71.27.-ha,75.10.-b,75.30.Ds 

o 
o 

Submitted to: New J. Phys. 



CN ■ 23 November 2011 

> 

m 



o 



% 



The time- dependent Gutzwiller theory for multi-hand Hubbard models 2 

1. Introduction 

The study of materials with medium to strong Coulomb-interaction effects has been 
a central subject for experimental and theoretical solid-state physics over many years. 
Despite enormous efforts and significant progress in some fields, however, our theoretical 
toolbox is still far from satisfactory for such systems. For quite some time, theoreticians 
in many-particle physics have focused on relatively simple model systems, such as the 
Heisenberg or the single-band Hubbard models. Only in the past ten years, attention 
shifted towards the study of more realistic models, e.g., multi-band Hubbard models. 
A very important impulse in that direction came from the limit of infinite spatial 
dimensions {D — >■ oo). The exact solution of Hubbard models in this limit leads to the 
Dynamical Mean Field Theory (DMFT), in which the original lattice model is mapped 
onto an effective single-impurity system that has to be solved numerically [H El |3l IH [5] . 
Although significant progress has been made in recent years in developing numerical 
techniques for the solution of the DMFT equations, it is still quite challenging and can 
be carried out only with limited accuracy. It is particularly difficult for the DMFT 
to study multi-orbital Hubbard models when the full (local) Coulomb and exchange 
interaction is included. 

An alternative method that also relies on infinite-D techniques is the Gutzwiller 
variational approach. It allows for the approximate study of ground-state properties and 
single-particle excitations with much less numerical effort than within DMFT and has 
been applied in a number of works in recent years [HI [3 El El [101 [HI [121 [131 [HI [IS [lEl [13 
[T8l [T9l [201 HH [221 [231 |2l]. Another approach that leads to the same energy functional 
for multi-band models is the slave-boson mean field theory [251 [SHI [271 [Ml [291 [30] . 

The theoretical interpretation of a number of experiments requires the study of 
two-particle response functions. For example, in magnetic neutron scattering the 
frequency- and momentum-resolved magnetic susceptibility is measured. The textbook 
method for the calculation of such response functions is the random-phase approximation 
(RPA), which can be interpreted as a time-dependent generalisation of the Hartree-Fock 
(HF) theory in the small amplitude limit, i.e., where the perturbation is considered 
to be sufficiently small. For electronic systems with medium or strong correlation 
effects, however, the ground-state description of a HF theory is well known to be often 
inaccurate. Therefore the RPA, as the time-dependent generalisation of the HF theory, 
is also a questionable approach for such systems. 

A time-dependent Gutzwiller theory for the calculation of two-particle response 
functions was developed for single-band Hubbard models by Seibold et al. [311 [32] . In 
recent years this approach has been applied with astonishing success to quite a number 
of such models and response functions [331 [Ml [351 [Ml [SH |38l [39l iQl SH [121 113] . It is the 
main purpose of the present work to generalise the time-dependent Gutzwiller theory 
for the investigation of multi-band models. A brief introduction into our method has 
already been given in Ref. [44j. All technical details, however, will be first presented 
here. 
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Our presentation is organised as follows: In chapters [2] and [3] we summarise the main 
results of the Gutzwiller variational theory for multi-band Hubbard models. In chapter H] 
the reader will be reminded of the derivation which introduces the RPA as a time- 
dependent generalisation of the HF theory. In a very similar way the time- dependent 
Gutzwiller theory ('Gutzwiller RPA') is introduced in chapterO The general Gutzwiller 
RPA equations are used in chapter [6] for the calculation of response functions for 
Hubbard-type lattice models. As a first application we study the magnetic susceptibility 
in a two-band model in chapter [71 A summary and conclusions close our presentation 
in chapter [HI The more technical parts of our derivation are referred to four appendices. 

2. Multi-band Hubbard models and Gutzwiller wave functions 

We study the general class of multi-band Hubbard models 

Here, the first term describes the hopping of electrons between N spin-orbital states 
0", a' on Ls lattice sites ?,j, respectively. The Hamiltonian 

fj ^ ^\ Sr^ 7y<Tl,<T2,<73,a4 2t gt A. A. I Y^ ^1,^2 At A. 
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Tl ,0-2, 0-3,0-4 



J,cri^*jO"2 



contains all local terms, i.e., the two-particle Coulomb interactions (~ Ui) and the 
orbital onsite-energies (~ ej). We further introduce the eigenstates |r)j of li\oc,i and the 
corresponding energies E'^\, i.e., 

-f^ioc,i|r)i = -Sr°i|r)i • (3) 

Within the Gutzwiller theory, the Hamiltonian ([1]) is investigated by means of the 
variational wave function 

|*G)=PG|^o)=n^«l^o)' (4) 

i 

where |\l/o) is a normalised single-particle product state and the local Gutzwiller 
correlator is defined as 

A = $^V,r'|r)..(r|. (5) 

r,r' 
For example, in case of the single-band Hubbard model 

«7^i o-=t,i i 

the local correlation operator reads 

Pi = K,d\d)ii{d\ + Ai,tlt)n(tl + \,U)^i{i\ + KiMi)i^{t\ 

+ hunhiii + KMum ■ (7) 

Here, we introduced the atomic states |r)j for doubly occupied sites \d)i, singly occupied 
sites |t)i and |J,)j, and for empty sites |0)j, as well as the abbreviation Aj^r for the diagonal 
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variational parameters Aj;r,r- In terms of the fermionic operators cfl, the operator ([7]) 
has the form 

where rijo- = cl^c^^. The correlation operator ([7]) is the most general Ansatz for single- 
band models without superconductivity. The latter would require additional terms of 
the form ~ |d)ii(0| and ~ |0)ii(rf|; c.f., Ref. [15]. 

3. Variational energy 

As shown in Refs. [TJ |16], the expectation value of the Hamiltonian ([1]) with respect 
to the variational wave-function (jlj) can be evaluated in the limit of infinite spatial 
dimensions. We consider the expectation values of the local Hamiltonian Hioc, and 
the one-particle Hamiltonian Hi separately in sections 13.11 and 13.21 The additional 
constraints, which arise through the derivation in infinite dimensions are discussed in 
section 13.31 In section 13.41 we recall how the standard Gutzwiller energy functional for 
a single-band model is recovered from our general multi-band results. 

3.1. Local energy 

The expectation value of the local Hamiltonian ([3]) in infinite dimensions reads [48] 

(ifioc)*G = / ^E^'^mr,r = E\oc , (9) 

r 

where 

mr,r' = ((|r)(r'|))^, = ( (pt|r)(r'|p))^^^ = ^ a* ^ a,, f,m?_f, (10) 



r,r' 



and 



m 





r,r' 



((|r)(r'|))^, . (11) 



To further evaluate the expectation value flTT|) . we introduce the basis of Fock states 
(i.e., 'Slater determinants') 

i/)=n^iio) (12) 

in which certain spin-orbit states a E I are occupied. Mathematically, the indices 
/ = (cTi, cr2, • • • , o"n) are considered as ordered sets of spin-orbit states a. Therefore we 
can use all standard set operations, such as / U a or I\a. In addition, we define the 
number of orbitals in / as |/|. The states |/) provide a basis of the local (atomic) Hilbert 
space. Hence, we can use them for an expansion of the eigenstates |r), 

\^) = Y.^i,r\I) (13) 
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and write the expectation value (ITT]) as 



m 



r,r' 



/,/' 



Finally, the uncorrelated expectation values of the transfer operators |/)(/'| 

ml,^{{\I){I'\))^,, 
can be written as the determinant 



m 



IJ' 



n 



jj' 



n 



J, J 



(14) 



(15) 



(16) 



Here, Qj i> are the matrices 



Qj 



( CO , 

01,0j 


•^01,0^ • 


• ■ ^oi,of^,| 


\ 


'^02,0; 


^^02,0^ • 


(^0 




\ ^°IIh-i 


0|i-|,02 


CO 

0|j|,of^, 


/ 



(17) 



in which the entries are the elements of the uncorrelated local density matrix 

CO^, = (c,Vc,,,)*„ (18) 

that belong to the configurations J = (ai, . . . , g\i\) and /' = (a^, . . . , crL,). The matrix 
0"^'"^ in (116!) is defined as 



^ 



J, J 



/ 1-CO 

■-^01,01 

02,01 


-CO 

01,02 

^ ~ ^02,02 ■ 




-CO \ 

01,0|j| 
02,0|j| 


V -C'o|,;|,01 


-CO 

0|j|,02 


. 1 


- C'o|j|,0|,,| / 



(19) 



withfJiG J=(1,...,A^)\(/U/'). 

In applications of the Gutzwiller theory to multi-band systems it would be quite 
cumbersome to evaluate the determinants (fT6|l if the local density matrix (TTSj) is non- 
diagonal. Fortunately, we are free to chose the local orbital basis in a way that suits 
us best. Therefore, we introduce an orbital basis, defined by (local) operators h^\ for 
which 



c 



7.7 



7.7' ~ "7,7'\"'i,7"'J,7')*o = '"'7 • 



With such a basis the expectation value (116!) has the simple form 



nil 



n„ 



(20) 
(21) 



76/ 7^/ 
Note that, for simplicity, we always use the same variable / for configuration states of 
the form ( 1T2|) irrespective of the underlying orbital basis (e.g., cj in ( [T2l) or h^ in ( 12T|) ). 
As will be shown in chapter |5l the time-dependent Gutzwiller theory requires to 
calculate the first and second derivatives of the energy with respect to all elements of the 
local density matrix, including the non-diagonal terms. Since the local density matrix 



The time- dependent Gutzwiller theory for multi-band Hubbard models 6 

enters the energy functional solely through matrices of the form (TTB]) we only need to 
expand these matrices with respect to small perturbations 

'y,j' ^ 7,7' "'" 7,7' l^^J 

up to second order in SC^^y around the diagonal ground-state matrix (pOl) . This 



expansion is explicitly carried out in Appendix A 



In our derivation of the ground-state energy all local onsite energies were considered 
as part of the local Hamiltonian ([2]). For later use, however, we also need an expression 
for the expectation value of a general local one-particle Hamiltonian 

^onsite = 5^ e'^'^'c^c,'. (23) 

This expectation value is given as 

a, a' 

where 

Q<.'= E ^krM,rA^Mc.'\r3)ml^r, (25) 

ri,r2,r3,r4 

is the 'correlated' local density matrix. 

3.2. Kinetic energy 

The expectation value of a hopping term in Hq, Eq. ([1]), is given as 



where we have introduced the (local) renormalisation matrix 

€'= E ^r.A^r3,r.(r2|al|r3)E^AA^;„r.</4- (27) 

ri,...,r4 hM 

The matrix HJ_^ j^ contains three different contributions depending on whether the index 
a' is an element of /in/4, /4\(/in/4), or J = (1, ... , iV)\(/iU/4). With the abbreviation 
fa,i = {I\clcJI) we can write Hf^j^ as 

Hjih = (1 - f.',iJ{h\dM U ^'X,/,u.' (28) 

+ {hV\c^,\Ii) [fa',hrn'j^\a',u + (1 - f'^'J4)^°/Xa',h) ■ 

The expectation value rrij^",, j in fl28|) has the same form as the one in (fT6|) . except that 
the index J has to be replaced by J\cr'. In case of a diagonal local density matrix one 
finds 

Ku = ShW,iAhV\cJh)Y^^^^ (29) 
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in agreement with results derived earlier [7]. Note that, in general, the renormalisation 
matrix is not Hermitian, i.e., it is 

c' ^ icr ■ (30) 

Using (l26l) . the expectation value of the one particle Hamiltonian Hq can be written 
as 

{Hoh, = L, J2 Ci(C^)*Ki,..,<,4 (31) 

where we introduced the tensor 

^.^.^i,<^^ = ;^E<r(Ci^^-^)*o- (32) 

3.3. Physical constraints 

As it turns out through the evaluation of expectation values in infinite dimensions, the 
variational parameters Ar,r' need to obey certain local constraints. These are 

{P^PU, = 1 , (33) 

{ciP^Pc^,)^, = (cic^,)^, = Cl,, . (34) 

Note that moving the operator P'^P relative to cj. or c^, in f p^ would not alter the 
whole set of constraints. With the explicit form (jS]) of the correlation operator P, the 
constraints read 

(35) 

r|ct|ri)x(r3|c,,|r>° . (36) 



3.4. Recovery of the 'standard' single-band energy functional 

In case of a single-band model, the atomic eigenstates |r) coincide with the configuration 
states |/). If we assume a translationally invariant ground state and the most general 
form of a local density matrix jlS] 

1^ (cic,)^„ {c[c,),, ) [ AO,, nl J ' ^'^^ 

where (Aj.^)* = A°^^ = A°, we find 

m% = {l-n',){l~nl)-\Ar, (38) 



1 


— 2-^ •^r,ri-^r,r2"^ri,r2 




r,ri,r2 


^0 


= Z^ -^ra.riVa.ra' 




r,r',ri,r2,r3 



m 



m 



= nl{l-nl) + \Ar, (39) 

°,. = K,. , (40) 

ml,=ny,-\Ar, (41) 

for those of the expectation values (llSp which are finite. Here we used the notation 

t=i and i=t . (42) 
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As a consequence, the expectation value of the local Coulomb interaction in (^ reads 
J2u{\d)u{d\)^^ = LM\Xd\'ml, . (43) 

i 

For the single-band model with the correlation operator ((Tj) and the local density 
matrix fl371) the elements of the renormalisation matrix have the form 



q: = A:A0(1 - nl) + A*A,4 + (A*A,,. + K,MK^ , (44) 

q: = Ai,{x:x, - x:x,) - a*a,,.< + x:^,x,{i - n°) . (45) 

Finally, the constraints in this case are given as 

1 = \X,\'m% + \X,\'m% + (|AtP + \X,M'H,, + (|A,p + \Xu\')ml^ 

+ {Xl,X, + X;X,,,)Al, + (A;,,At + AtA,,t) a;,, , (46) 

nl =i\X,\' + \X^M'X4+\X,\'ml., (47) 

A0_, = - (A;^A. + A;A^,.)m°,, + |A0pA°_, . (48) 

As mentioned before, it is possible to overcome the complications that arise from a 
non-diagonal local density matrix by a simple transformation 

h\ = J2^^A (49) 



to a new orbital basis for which the local density matrix is diagonal by definition, 
Chlh,)^, {h\h,)^^ \^(n^, \ 

In this new basis the constraints have the rather simple form 

1 =Xlm% + XJml, + Xlml^ + Xlm%, (51) 



^0 



n, 



., A>;,, + Xlm% . (52) 

where the non-diagonal constraints are automatically fulfilled by working with a diagonal 
correlation operator (Ai_2 = 0). Equations (1^ - (132]) can be readily solved by introducing 
the expectation values 

fhd = ~Xlmla , (53) 

m0 = A0m0 = l-hl-h^ + rhd , (54) 

m^ = A^m°_^ = n?^-md , (55) 

which leaves us with only one variational parameter, the expectation value m^^ for 
a double occupancy. The resulting renormalisation matrix is then diagonal and its 
elements have the well known form Il9l [501 El] 



q-y = q'^ = A0A^(l-n°)+AtjA^n° = (^/^f^^ + ^/m^^^ .(56) 



nO(l-nO) 



Hence, the single-particle energy ( 13T|) is given as 

7 «7^i 
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where we used the orthonormahty relation 

y^ Ua,-yK,-y' = ^1,Y ■ (58) 

(T 

In order to formally show the equivalence of both approaches we write the 
parameters Ar,r' in © as 

A0 = A0 , (59) 

Xd =Xd, (60) 

K,ct' = ^Ma,7<',^\ • (61) 

7 

With these relations it is easy to show that the constraints (l46 l) -f H8l) are indeed fulfilled. 
For example, the first constraint (H6|) can be written as 

1 = Xlm% + Xlml, + J2 K",.K", .'<,.' ■ (62) 

If we use 

^l,a' = ^Kr,^'^'n'^^^, (63) 

7 

and the orthonormahty relation (I58p we readily find that 0621) is indeed solved by the 
parameters (l59l) - fl6Tl) . In the same way one can show the equivalence of the ground-state 
energy functionals. 

4. The Time-Dependent Hartree-Fock Approximation 

The approximation most frequently applied to two-particle Green's functions is the 
'random-phase approximation' (RPA). This approach can be derived in various ways, 
e.g., by an equation of motion technique or in diagrammatic perturbation theory [52] • In 
this section, we use a different derivation which introduces the RPA as a time-dependent 
generalisation of the Hartree-Fock theory; see, e.g., Refs. [S31 HI]. If derived in this way, 
the approach can be generalised quite naturally in order to formulate a time-dependent 
Gutzwiller theory. This will be the subject of chapter O 

4-1. The Hartree-Fock approximation 

In the Hartree-Fock approximation a single-particle product wave function |\E'o) is used 
in order to investigate the ground-state properties of a many-particle system. Note 
that such wave functions are included in the Gutzwiller variational space by setting 
-^j;r,r' = ^v,v'- The expectation value of a many-particle Hamiltonian with respect to 
a Hartree-Fock wave function is a function of the single-particle density matrix. For 
example, for the Hamiltonian ([T]) it reads 

E«^(p) ^ {H)^, ^ (64) 
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where 

PUa'Uia) = (4<xS,a')*0 (65) 

are the elements of the single-particle density matrix p and 

^lo^i(p) = 2 Z2 Ui^''"^'"^'"^ P{ia4),(i<7i)P(i<73),(i<72)-P{ia3),{mi)P{m4),{ia2) (66) 

CT1,CT2, 

CT3,o-4 

is the expectation value of the Coulomb interaction in the Hamiltonian ([2]). Note that it 
will be more convenient both in the time-dependent Hartree-Fock and Gutzwiller theory, 
to use a different order of subscripts in the definition (1651) of density matrices than, e.g., 
in section |3] or in previous work on the Gutzwiller theory. 

To keep notations simple, we use the abbreviations v = {i, a) for local single-particle 
states and Y = {v, v') for pairs of these indices. For example, the elements of p can 
then be written as 

PY = Pvi,V2 = P(Jl,ffl),(J2,ff2) • (67) 

With these new notations, the Hartree-Fock energy (!64l) reads 



-^ \P) / J ^Vl,V2Pv2,Vl ~r r) / J Pv4,Vi'^{vi,V4),{v3,V2)Pv3,V2 

= Yl ^^P^ ^ 2 ^ PyWy,y'Py' (68) 



Y,Y' 



where 



and 



^(m),(i'T2) = C-'"^' + f^ijer'"' (69) 



TI/ — 7'ra"i,o-2,a"3,o-4 rro-i ,0-2,0-4,0-3 /yn\ 

\^{vi,V4),{vs,V2) = l^i - l^i UUj 

for indices Vk = {i, <Jk) that belong to the same lattice site i. Further, we introduced 
the 'inverse' index Y = {v',v) for Y = {v,v'). Note the symmetries 

'^^{v-i,V4),{v3,V2) '^ {V2,vs),{v4,vi) '^^ {vi ,V3) ,{V4 ,V2) 5 ['^) 

which will be employed in the following section. 

The energy functional (168|) has to be minimised with respect to all density matrices 
which belong to a single-particle product state. Such matrices are idempotent, i.e., they 
obey the matrix equation 

p' = P- (72) 

If one imposes this constraint via a Lagrange parameter matrix fj with elements t]^^^> , 
the following equation has to be solved 
d 



dpv. 
This condition leads to 



E^'Cp)-tTif,ip^-~p) =0. (73) 



h{p) -\- fj — TIP — pfj = (74) 
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where we introduced the matrix h{p) with the elements 

hy(p) = /t^^^Ip) =^y + Yl ^y,y'Py' • (75) 

Equation frM|) is solved if p satisfies both fl72l) and 

[h{p),p]=0. (76) 

Starting with a certain density matrix p we can introduce the 'Hartree-Fock' basis 

\a) = ^u^,a \v) (77) 

V 

of states which diagonahse the Hamilton matrix h{p), i.e., 

y^ hv,v'{p)Uv',a = EaUy^a ■ (78) 



Equation (!76|) is then solved by setting 

Pa,a' = 5a,a'Q{EF " ^a) (79) 

where the Fermi energy Ep is determined by the total number of particles 

N = J2^iEF-E^). (80) 



The density matrix (1791) has to be reinserted into fl68l) . fl75|) until self-consistency is 
reached. We denote the solution of these equations as p^ and introduce the corresponding 
Hamilton matrix 

/i° = h{p'') . (81) 

4-2. Equation of Motion for the Density Matrix 

We consider two-particle Green's functions of the form 

G'(.„«,),K,«4)(i - t') = ((4,(i)c„,(t); tl^it')c^,it'))) (82) 

^ -ie(t-t')($o|[al,(t)c„,(t),ct3(t')c„,(t')]|$o) , 

where |$o) is the exact ground state of our multi-band Hubbard Hamiltonian ([1]), and 
Cv (t) is the Heisenberg representation of the operators ci, with respect to H. As shown 
in most textbooks on many-particle physics, the Green's functions (!82l) naturally arise 
in 'linear-response theory' because they describe the time-dependent changes 

/oo 
^t'G{v2,v^),{vs,V4,){t - t')fvi,vSt') (83) 

of the density matrix p in the presence of a small time-dependent perturbation 

Vfit) = J2U'it)^l^' (84) 



v,v' 
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added to H [SHI \5E\ E7] . After a Fourier transformation and using again the abbreviation 
Y = {v,v'), Eq. (183]) reads 

5pY{uj) = J2GY,Y'i^)fY'iuj) (85) 

Y' 

with 



oo 



Gy,y'{^)= drGy,y/(r)e^^", (86) 



oo 



and fri^) and 6py{uj) defined accordingly. 

Ideally, we would like to calculate the time dependence of the density matrix 

p^,,^{t)^{^{t)\cicj^{t)), (87) 

where |^(t)) is the exact solution of the time-dependent Schrodinger equation for the 
Hamiltonian 

H{t) = H + Vf{t) . (88) 

The expectation value flHTl) obeys the Heisenberg equation 

-ip.,,„(t) = (vl/(t)|[^,ctcj|vl/(t)), (89) 

which contains the commutator 

[Hit), clc^,] = J2 i^v„. + fvU^MA' - E (^-'.-1 + f-''-^ i^M^i (90) 

VI VI 

' n / J y^ {v\,v-i),{v,V2)'^vx^v-i^v''^v-i^ ^ {v\,V2).,{v-i,v')'^VY^v^v-i^vz] ' 

In the time- dependent Hartree-Fock approximation, it is assumed that the solution 
|\I'(t)) of the Schrodinger equation at any time t is approximately given by a single- 
particle product wave function. In this case, the expectation value of the commutator 
(!90|) can be evaluated by means of Wick's theorem. This leads to the equation of motion 

im = [Hp{t)) + f{t),m] (91) 

for pit), where the matrix h{p) has been introduced in (1751) . Equations (1751) and (I9T1) will 
be crucial also for our formulation of a time-dependent Gutzwiller theory in chapter [S] 

4 ■ 3. Expansion for Weak Perturbations 
We are only interested in cases where 

Vfit) ^ 5Vjit) = 5^5/.,.,(t)ctc„, (92) 

v,v' 

is a weak perturbation to the time-independent Hamiltonian H. In this case, the density 
matrix p{t) and the Hamilton matrix h{t) are given as 

Pit) ^ p° + 6pit) , (93) 

hit) ^h^ + 6hit) , (94) 
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where Sp(t) describes a 'small' time-dependent perturbation around the ground-state 
density matrix p^, and 

hy =^y + Yl ^y,y'Py' , (95) 

Y' 

Shrit) = ^WY,Y'SpY'{t) . (96) 

Y' 

With the expansion (!93|) -( 1M|) . the equation of motion ( 19T|) becomes 

=[/^°,pr (97) 

i5p{t) = [h', 5p{t)] + mt) + 5f{t),p'] . (98) 

These equations have to be solved for density matrices p{t) that obey the matrix 
equation (1721) . After applying the expansion ( l93ll . Eq. (1721) reads (to leading order 
in Spit)) 

p' = {P'T . (99) 

6p{t) = p'dpit) + 5p{t)p'' . (100) 

Note that Eqs. (1971) . (I99l) just recover the time- independent Hartree-Fock equations 
derived in section 14.11 

4.4- Random -phase approximation (RPA) equations 

Mathematically, the density matrix is a projector onto 'hole'-states, ph = P°. In addition, 
we define the projector onto 'particle'-states as 

Pp = 1-P°. (101) 

With these two operators, we can decompose all matrices into their four components 

dp'^^it) = pJpit)Pn. , (102) 

6r'"{t)^pjf{t)p^, (103) 

h'-^^^ = p.,h^p^ , (104) 

where v,w ^ {P;h}. Note that h^'^^^ has the elements 

/i°;^T = S^^y,6a,a'Eo, . (105) 

An evaluation of the condition (llOOp for the components 5p"^{t) yields 

5p^'"(t) = 5p^^(t) + 5p™(t)5p^'"(t) + 5p'"^(t)5p'""'(t) , (106) 

and 

,5p«'"'(t) = (107) 

where v ^ w. Hence the components 5p^'^{t) and 5p^^{t) can be neglected in the 
following compared to the leading fluctuations 5p^^{t) and dp^'^it). 

We express the time-dependent quantities 5p"^{t) and 6f^^{t) by their respective 
Fourier transforms 6p^^{uj) and 5f'"^{uj). The equation of motion fl98|) then leads to 

+ u5pll^^{uj) = {E^, - E^,)6pll^,^{u) ± (Shll^^ioj) + 5fZ^^{u)) (108) 
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where the plus and minus signs correspond to vw = ph and vw = hp, respectively. With 
the abbreviation A = (ai, 02) for pairs of indices a we find 

6hT{iu) = -J2 Ua,a'{^p7{^) + ^p7{^)) ■ (109) 



Here, the elements of the matrix U are given as 



Ua,. 



U, 



(ai,a2),(«i,a2) 



Yl '^*vi,ai'^V2,a2^(^u^2),K,v'^)Uy'^,a[ 



Uvi„ai, 



.(110) 



"[•"2 



The coefficients m„ q in flllOp have been introduced in Eq. fl77|) and determine the 
solutions |a) of the Hartree-Fock equations. Equations (11081) and (11091) then yield 



{u-E) 











+ u 






\5p 



with a matrix E defined as 



Ea,A' — -S(ai,a2),(a'i,a^) — ^ a^ ,a[^ a2 ,a'^{E a^ — Ea^) . 



By comparing Eqs. ( Illip and ( |85l) we find 

for the inverse of the two-particle Green's function 



u + i6-E)\ J _°^ \+U 



(111) 



(112) 



(113) 



Ga,A'{(^) = G(^aua2),{a[,a'2){(^) 



E 



'^v\,a\'^V2,a2 



(i>i,i;2),K,»^^)l^)«t,i,a'l"f^,a^ 



(114) 



V-i_,V2,V\,V'^ 



Here we have added an increment i5 with 5 = 0^ in order to ensure the correct boundary 
conditions of a retarded Green's function. For t/ = 0, the inverse Green's function 01130 
reads 

V-\u) = ±{uj + i5-E) (115) 

which leads to 



^A,A'{,^) = ^{ai,a2),(a'^,a'^){^) — ^ai,a'i'^Q:2,c 



r 02, 0:2 






'■uj-{E^, -Ea^) + i6 ■ 



;ii6) 



Note that F is not the exact Green's function for the single-particle Hamiltonian Hq since 
we just set f/ = in (11131) . but kept finite the 'Hartree-Fock self-energy' contributions 

J:a = J2Wa,a'p''a' (117) 

A' 

which usually change the 'eigenvalues' Ea in (I116p : c.f. Eqs. fl75]) and f l75]) . 
With the Green's function (I116P we can write (I113P as 

G{u) = r{u)[i + ur{u)]-^ (118) 

= t{io) + t{io)UG{io) (119) 

where, in the second line, we expanded [1 + Ur{u)]~^ into a power series with respect to 
f/F. Both Eqs. (I118p . (lll9p are familiar expressions for the two-particle Green's function 
in the random-phase approximation. 
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5. Time-Dependent Gutzwiller Theory 

The time-dependent Gutzwiller approximation has been first introduced for single-band 
Hubbard models by Seibold et al. [311 132] • In this section, we generalise this approach 
for the investigation of multi-band models. To this end, we set up an effective energy 
functional of the density matrix in section [5?T| which is used in sections l5 .2115 .41 to derive 
the Gutzwiller RPA equations. 

5.1. Effective Energy Functional 

As summarised in chapter [3l the expectation value of the multi-band Hamiltonian ([1]) 
in the Gutzwiller theory is a function of the variational parameters Ar,r' and of the 
one-particle wave function |\E'o)- Like in the Hartree-Fock theory, the single-particle 
wave function |\E'o) enters the energy functional solely through the elements fl65l) of the 
non-interacting density matrix p. It is therefore possible to consider the energy 

E = E{\,p) (120) 

as a function of the density matrix p and of the 'vector' 

A = ({A^,r'},{Ar,r'}) = (Ai,...,A„J (121) 

of Up variational parameters Ar,r' (and Afp, for F ^ F'). The density matrix in 
the energy functional (11201) must be derived from a single-particle wave function and, 
therefore, it has to obey the condition (1721) . Note that, in the following considerations, 
the density matrix will either be considered as a matrix (with respect to its two indices 
(i, 0") and (j, cr')) or as a vector (with respect to its single index Y). To distinguish both 
cases, we will denote the density matrix p in some equations as p in order to indicate 
its vector interpretation. 

The constraints (I35l) -( l36l) are also functions of A and p and will be denoted as 

(?„(A,p)=0 , l<n<n,. (122) 

Here, ric is the (maximum) number of independent constraints, which, due to 
symmetries, is usually smaller than its maximum value iVg^ + l, where Nso is the number 
of spin-orbital states per lattice site. We assume that the functions (I122p are real, i.e., in 
case of complex equations (I33l)-(l34l) their real and imaginary parts are treated separately. 
By solving Eqs. (11221) we can, at least in principle, express n,, of the variational 
parameters (= A^) through the density matrix py and the remaining 'independent' 
parameters (= A'^), 

Ai = Ai(A\p). (123) 

In this way, we obtain an energy functional 

E«^(A',p-) ^ E{X^{X\p),X\p) . (124) 
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which has to be minimised without constraints apart from Eq. f l72p and the condition 
that the total particle number 

N = J2pv,v (125) 

V 

is conserved. 

For a fixed density matrix p, the minimisation of f ll24p with respect to the 
parameters A^^, 

•9 -nGAMi 



^-E^^iX\p)=0, (126) 

determines these parameters 

y = m (127) 

as a function of p. This allows us to define the 'effective' energy functional 

E^^ip)=E^\X\p),p) (128) 

which, for a fixed density matrix p, is given as the minimum of E^^ with respect to A'. 
With this effective functional we will formulate the time-dependent Gutzwiller theory 
in the following section. 

Using a Lagrange-parameter matrix fi as in chapter 14.11 we find 
d 
dpv. 
which leads to 

0=[/i(p),p]. (130) 



E^\p)-ii{f,{~p'-p) 



= (129) 



Here we introduced the matrix h{p) with the elements 

nr(,) - ^ . (131) 

OpY 

and used again the notation Y = {j,a']i,a) for Y = (z, a; j, cr'). The self-consistent 
solution of Eqs. fll30p - fll31|) then yields the ground-state density matrix p°, the matrix 
/i° = h{p^), and the corresponding single-particle 'Gutzwiller-Hamiltonian' 

5.2. Gutzwiller RPA Equations 

The derivation of RPA-type equations within the time- dependent Gutzwiller theory 
goes along the same lines as discussed in chapter H] for the time-dependent Hartree-Fock 
theory. We add a small time-dependent field 

5Vf{t) = Y, SfLj,.'it)i.c^,.' + h-c. (133) 

i,jr';(T,cr' 

to our multi-band Hamiltonian ([T]). With the particular time dependence 

Sfl.,y{t)=6fl.^^y{u)e-'-' (134) 



The time- dependent Gutzwiller theory for multi-hand Hubbard models 17 

the expectation value of 6V{t) reads 

Ef{p)= Yl <^A-;.v(^)e~'"V,V;.,<x + c.c. , (135) 



i,j;a,a' 

where 



5f\^,,,^,{u) = S,Jfl,,Uu^)^^^^^^ (136) 

Pi,<J2;i,cri 



+ (i-ME^^^i..4H<ite 



r[,ai. 



The renormahsation matrix q and the (correlated) local density matrix C*^ are defined 
in equations ( 127|) and fl25l) . respectively. With Eq. (I127P they can both be considered as 
functions of p. 

The time-dependent field induces small fluctuations of the density matrix, 

Py = Py + Sprit). (137) 

Our main assumption is now that Spvit) obeys the same equation of motion, 

i6m = [/i°, 5pit)] + mt) + 5/(t), pO] , (138) 

as the density matrix in the time-dependent Hartree-Fock theory; see Eq. ( p8|) . Here, 
however, the Hamilton matrix 

h{t) ^ h%t) + 5h{t) (139) 

is not derived from the Hartree-Fock functional fl68|) . but from the effective energy 
functional f lT^ . 

hvit) = ^^^^(p) ^ 4 + y KYY'hY'{t) = h''y + 6hY{t) , (140) 

dpY V ' 

where the matrix K is given as 



K 



¥,¥' 



OpyOpY' p=pO 

The diagonalisation of h^ (or equivalently of the Gutzwiller Hamiltonian h^) yields 
a basis \a) with 

hlc' = h\ = <a'E^ (142) 

and a ground-state density matrix that is given as 

Pla' =Pa = Sa,a'Q{EF - E^) . (143) 

With the projectors ph = p° and pp = 1 — p°, we define the particle and hole components 
of all matrices, as we did in Eqs. (I102p - (11041) . The components 5p"^{t) of the density- 
matrix fiuctuations obey Eqs. (I106p -( TT07|) . i.e., to leading order we can neglect 5p^^{t) 
and 6fF^{t). Hence, after a Fourier transformation we end up with the same form of 
RPA equations. 



(^-^)| I ', \+K 



( 5p^\u) \ / 5~P\l 



{6p^^{uj)) {6Mu;) • ^^^^^ 
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as in Eq. fillip . Here, however, the bare matrix of Coulomb parameters U is replaced 
by the matrix K, defined in fll4ip . and the energies E^ in the matrix E, Eq. f lll2p . are 
the eigenvalues of the Gutzwiller Hamiltonian f ll32p . The comparison with (185|) leads 
to the final result 



G{uj) 



(w + i(5 - E) I ; ^^ \+K 



(145) 



for the two-particle Green's function matrix within the time- dependent Gutzwiller 
approximation. 

One should keep in mind that the external 'fields' Sf{uj) in fll44p are 'renormalised', 
i.e., they are not the bare fields as they appear in ( 185!) . see Eq. ( 11360 . On the other 
hand, on the l.h.s. of ( l85l) appears the 'correlated' expectation value of the density 
matrix, while in f ll44p we work with the fiuctuations of the uncorrelated density matrix. 
Therefore, the 'true' Green's function seen in experiments may, in fact, be given as 

Gy^Y'i^) = cy,y'Gy,y'{uj) , (146) 

with certain frequency independent factors Cy^y- These factors, however, are of minor 
importance since they only affect the overall spectral weight and not the frequency 
dependence of the Green's function matrix G{u)). We can calculate them with the 
assumption that correlated and uncorrelated density-matrix fluctuations are related 
through the same renormalisation factors as the corresponding ground-state density 
matrices. For the one-band model it has been checked that this prescription is in fact 
the correct procedure for which the correlation functions fulfil the standard sum rules 
[3111371113]. 

5.3. Second- Order Expansion of the Energy Functional 

For an evaluation of the Gutzwiller RPA equations (I144p . we need to determine 
the matrix K which is given by the second derivatives (I14ip of the effective energy 
functional (I128p . To this end, we expand E^^ up to second order around the ground 
state values p° and A''° = A'(p°), 

E'^^y, p) = Eo + tiCh'Sp) + ^[Y.6pYM^y,6pY' + J2 ^^'zM^%'^X'z' 

¥,¥' Z,Z' 

+ J2 (S>}zMz'^y6py + 6pyM^^z^X'z 
Z,Y 

= Eo + iiih^dp) + (5^(2) . (147) 

Here, we introduced the matrices M^p, M^^, M^^, M^'^ with the elements 

OpyOpyi 
a2E;GA 
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where the second derivatives on the r.h.s. are evaluated for p = p^ and A' = A''°. Note 
that there is no hnear term ~ A'^ in fll47p because of the minimisation condition fll26p . 
For our further evaluation, it is useful to write the second order terms in Eq. ( I147p in a 
more compact form by means of matrix-vector products, 

1 



^^'" 2 



{6p)^MPP6p + 2{6XYM^P6p + {6XYM^^6X' . (151) 

Here we used the symmetry 

M^p = [MP^f . (152) 

In the effective energy functional (I128p the parameters A^ are determined by 
the minimisation condition (I126p . Applied to our second-order expansion (I15ip this 
condition yields 
d 



dSXz 



-SE^^\5X\Sp) = 0, (153) 



which gives us the multiplet-amplitudes 



5X' 



M 



x\ 



-1 _ 



M^P6p (154) 



as a linear function of the densities 6p. This result leads to the quadratic expansion 

E^'^(p° + 5p} = Eo + iiih^p) + ]^{5p)^k5p, (155) 



k = MPP - MP^ 



M^^ 



-1 _ 



M^P , (156) 



of the effective energy as a function of the density fluctuations 5p. In earlier work on 
the time-dependent Gutzwiller theory, Eqs. (I153p and (I154p have been denoted as the 
'antiadiabaticity assumption'. In fact, these equations have the physical meaning that 
the local multiplet dynamics, described by fluctuations 5A^2(t), are fast compared to 
those of the density- matrix fluctuations 5py(t). We will use the phrase 'antiadiabaticity 
assumption' in this work too although, strictly speaking, in our derivation it does not 
constitute an additional approximation. 

With the functional (11550 . we could now proceed with our evaluation of the 
Gutzwiller RPA Eqs. (11440 . For practical applications, however, it is more convenient 
to determine the 'interaction kernel' (I156P in a way that avoids the explicit solution of 
the constraint equations (I122p . This alternative procedure is the subject of the following 
section. 

5.4- Lagrange-functional expansion 



In the second-order expansion, described in section 15. 3[ we implemented the 
constraints (I122p by explicitly eliminating a certain set of Uc variational parameters. 
Although such a procedure can, at least in principle, always be applied, for the numerical 
implementation it is more convenient to impose the constraints by means of Lagrange 
parameters. To this end, we define the 'Lagrange functional' 

ric 

L(A, p. A) = E{X, p) + Y, ^n9n{\ P) (157) 

n=l 
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which depends on all variational parameters A, the density matrix p{=p) and the n^ 
Lagrange parameters A„. The optimum variational parameters A^, density-matrix 
elements py-, and Lagrange parameters A° are then determined by the equations 



dL 



dX, 



dL 



OK 



dL 



dpy 







(158) 



A=AO,A=AO,p=pO 

which have to be solved simultaneously. 

We expand the Lagrange functional to leading order with respect to parameter 
{6\z, SAn) and density fluctuations (Spy)- The second-order contribution has the form 

5L(2) = 1 ^ 6pYL'^'y,6py, + J2 ^^ZL^Y^PY + \Y. ^^zLf^z'^Xz' 



Y,Y' 



zx 



Z,Z' 



+ 5^mJ5^^5Az 



dXz 



^t^^ 



(159) 



Z ^ Y 

with matrices L^^, L^p, L^^ defined as in Eqs. (I148I) -( TT501) only with E'-^^ replaced by L. 
The antiadiabaticity conditions 
d 



d5\z 
d 



dSAn 
yield the Uc equations 



5L(2) = , 
5L(2) = , 






(160) 
(161) 

(162) 



and the Up equations 



dgn 



J2 L'z'z'S^z' + E L'/^y6pY + Yl ^^^n = . (163) 

Z' Y n ^ 

Together these equations allow us to express the n^ + n^ parameter fluctuations (5A„, 
SXz in terms of the density fluctuations 6pY- These can be reinserted into (11591) to 
obtain the desired quadratic functional solely of the density fluctuations, 



5-Z^ = - y ^5pyKy,y'^Py' 



(164) 



Y,Y' 



In Appendix B.l, we prove that the interaction matrix Ky,y' iii (I164p is, in fact, identical 
to Ky,y' in Eqs. (1T55D-(TT56D. 



6. Two particle response functions for lattice models 

In the previous chapter we have developed the general formalism of the time-dependent 
Gutzwiller theory for the calculation of two-particle Green's functions. We will be more 
specific in this section and explain in detail how the response functions which are of 
interest in solid-state physics can be calculated within our approach. 
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6.1. Two-particle response functions 

In solid-state physics one is usually not interested in the full two-particle Greens- 
function G as it has been defined in ( l82l) . The properties, relevant for experiments, are 
certain linear combinations of elements of G. For our translationally invariant model 
Hamiltonians ([I]) these are in particular the two-particle response functions 

G(a,,.,u<rs,..)iR^-RJ,t-t') ^ ( (c,V, (t)c,,,, (t) ; c],,3 (t')c^. ., (f) )) (165) 
or, more importantly, their Fourier transforms 

oo 

— oo '-^ 

^ "^ ^\\'^fc,aiCfc+5,o-2' %+<J,o-3%,o-4)/^ ■ (166) 

^ k,k' 

Here, we introduced the fermionic operators 



A(t) _ 



1 ^,Ti«V.^£(t) (167) 



and the usual notation 

((6; 6'))^= / rfr((O(r);O'(0)))e^'^- (168) 



for the Fourier transform of a Green's function with arbitrary operators 0,0' . With the 
abbreviation v = {a, a') for spin-orbit indices and the operators 

^l^^L^^^4rJ2^UA+,,.. (169) 

we can write (I166p as 

G.y{q,u) = {{Al;{Al,y))^. (170) 

The Green's functions (I166p are still quite general since they include all possible 
channels of local coupling o"i ^ o"2, o"3 ^ c^^. In experiments one usually measures 
response functions which are certain linear combinations, 

Ge{q,uj) = y^^KyGv^^'{q,U})Ky> (171) 

v,v' 

of some of the Green's functions (I166p . defined by the matrix k,^ = n^y- For example, 
the transversal spin-susceptibility x(^) ^) is given as 

x{q,u) = ^{{S^;SZ^))^ (172) 



-'s 



where 



i k b 

J k b 
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are the usual spin-flip operators and b is an index for the orbitals at each lattice site i. 
The spin susceptibility of a two-band Hubbard model will be investigated in chapter [3 

6.2. Response functions in the time- dependent Gutzwiller approximation 

In order to apply the time-dependent Gutzwiller approximation, as developed in 
chapter [5l we have to expand the Lagrange functional (11571) up to second order with 
respect to density-matrix {6p) and variational-parameter fluctuations (5Ar,r')- This 
means that we need an expansion of the constraints fl35l) - fl36l) . of the local energies iQ- 
( ITTj) and (El]), and of the kinetic energy (l3T1) -( l32l) . The second-order expansion of the 
kinetic energy is more involved than that of the local energies and of the constraints. 
In the latter there are only contributions from fluctuations at same lattice sites while 
in the kinetic energy local and non-local fluctuations (such as <^(c|o.c,o-')*o) couple. 
Nevertheless, the calculation of the second-order Lagrange functional is tedious but 



otherwise straightforward. We therefore refer to Appendix C where the details of this 
derivation are presented. As shown in that Appendix, it is useful to introduce the 
operators 

1 



^w — ^a,,CTo,a'a' " /7^ Z^ ^fc ' ^ '^k,a'Sk+n,ai ' i^'^^) 



■'s 



^w — -Dcri,(T2,ai,4 — n^Z-^'^k+q ^k,a'^'^k+q,a[ ' (1' ') 

and to define the auxiliary Green's function matrix Il{q,uj) with the elements 

f {{Al;{Al,n. {{Al-{Bl,n^ {{Al-{Bm^ 



liv ,v' (q, w) 



{{Bl,{Kn. {{Bl-{Bm^ {{BI;{bI)% 



\{{K;{Knu. {{B^-{Bl,)^)% {{B^;{Bj))J 

(178) 

We are actually interested only in the first 'element' of this matrix, i.e., the Green's 
functions (11701) since they allow us to determine any response function of the form (I17ip . 



As shown in Appendix D , however, the time-dependent Gutzwiller approximation leads 
to the following equation for the entire matrix (I178p from which (I170p can be extracted, 

II(g, uj) = {l + II°(g, uj)V'^)-'ti\q, to) . (179) 

Here, V^ is the effective second-order interaction matrix, introduced in ( JC.401) . and 
H°(g, u) is the Green's function matrix (I178P evaluated for the single-particle Gutzwiller 
Hamiltonian (I132p . As shown in Refs. [TOl HB], this Gutzwiller Hamiltonian h^ = H^^ 
for our lattice Hamiltonian ([T]) has the form 

^o'' = E E (^T"^ + v..,..)ci.A^, = E E Ef^^^KAa (180) 

k cri,(T2 k a 

where the Lagrange parameters ?7o-i,o-2 ^"^^ determined by the minimisation of the 
variational ground-state energy and ej}''^' is defined as 



^k 



= ECi(Cp*<'^"^- (181) 
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The creation and annihilation operators /i^^, of the effective single-particle 
Hamiltonian f llSOp can be written as 

^L=E<X<.' (182) 

ka = J2K^y%^ ' (183) 

where the coefficients Uo-,« are determined by a diagonalisation of f llSOp . With these 
eigenstates the calculation of Il°(g, cj) is now a simple task. For example, the first 
element ((A^; (A^,) ))° is given as 

{{^lu^.-^iK'J)Z (184) 



T Z-^ / . \\'^k,a9.'^k+p.cn'j '''k'+a.a'/''k'.au)^ 



)^ iv^ Yv^^'^ iv^'^'^ Yi/- 

^ k,k' °i'"2 



k «i,02 



T!l \ a-^,OLi' a^ 



k 

(72,02'' "'cri,ai\"'a'^^,ai' "'0-2,02/ Q 



00 - (E,+,,,, - i^fe,. J + i<5 ^""''^^ """^'^^^ ' 



where 



nl^ = eiEp - Ek,a) (185) 

is the ground-state distribution function (11431) . In the same way, we can calculate all 
other elements of 11*^ (g, u) . The result is always the same as in fll84p only with additional 
factors ~ e^'"^ or ~ e^,'^ due to the definition of the operators fll76l) - fll77p . For example, 
the second element in fll78p leads to 

{{^luo2-ABlo..',,.'M (186) 

_ ^ Sr\^ ^"^2,02^ "oi,Oll"a^,aiJ V."2 03,04/ _n X 

T 2-^2-^ ,, ( TP TP "i _L iA '^ l'^fc,a2 "'k+q,ai) ■ 

k OL\,Oi2 • H, L , Al 

To summarise, with Eqs. fll79p . f ll84p . f ll86p . and the interaction matrix (IC.40p we 
are now in the position to investigate any two-particle response function for our general 
class of multi-band models ([1]). As a first example, we study the magnetic susceptibility 
for a two-band model in the following section. 

7. Magnetic susceptibility of a two-band Hubbard model 

In this chapter we investigate the magnetic susceptibility of a two-band Hubbard model 
in three spatial dimensions. The model Hamiltonian and the Gutzwiller wave functions 
which we use for its investigation are introduced in section mi In section [72] we discuss 
the Green's function matrices that we need to study in order to calculate the magnetic 
susceptibilities within the RPA and the Gutzwiller-RPA schemes. The numerical results 
for the two-band model are presented in section 17.31 



The time- dependent Gutzwiller theory for multi-hand Hubbard models 



24 



# 


Atomic eigenstate F) 


Symmetry 


energy Er 


1 


lt,t) 


'A, 


U' -J 


2 


(lt,i)+ ;,t))/v^ 


'A, 


U' -J 


3 


114) 


'A, 


U'-J 


4 


(lt,i)- iA))/V2 


'E 


U' + J 


5 


(ti,o)- o,n))/V2 


'E 


U-J 


6 


(|ti,0)+ 0,ti))/v^ 


'A, 


U + J 



Table 1. Two-particle eigenstates with symmetry specifications and energies. 



7.1. Model and variational ground state 

We investigate a Hubbard model with two degenerate Cg orbitals per site on a cubic 
lattice. The local Hamiltonian (jj]) for this system can be written as 

H^^ = U^ ne,tne,i + U'^ hi^sh2,s' - J^^ ni^sn2,s (187) 

e s,s' s 

s 

Here, e = 1, 2 labels the Cg orbitals, s =t, i is the spin index and we use the convention 
f =J,, \ =-\. Due to the cubic symmetry the Coulomb parameters f/, U' and the 
exchange parameter J are related to each other through 

U' = U-2J. (188) 

Hence, only two of these three parameters can be chosen independently. 

There are four spin-orbital states a = (e, s) per atom, leading to a 2^ = 16- 
dimensional atomic Hilbert space. All eigenstates |r) of H^^ with particle numbers 
N y^ 2 are simple Slater determinants of spin-orbital states \a) and their energies are 

Er = {N = 0,1) , 



Er = U + 2U' -J {N = 3) , (189) 

Er = 2U + AU' -2J {N = A) . 

The two-particle eigenstates are slightly more complicated because some of them 
are linear combinations of Slater determinants. We introduce the basis 



|s,s) 

lti,o) 

|0,ti) 



cl,sct' |0) 
cl,/i,;|0) 



|0) 



(190) 
(191) 
(192) 



of two-particle states, which are used to set up the eigenstates of H\oc;i, see table 17. ll The 
states of lowest energy are the three triplet states with spin S = 1, which belong to the 
representation A2 of the cubic point-symmetry group. Finding a high-spin ground state 
is a simple consequence of Hund's first rule. Higher in energy are the two degenerate 
singlet states of symmetry E and the non-degenerate singlet state of symmetry Ai. 
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Figure 1. Left: Model density of states at the Fermi energy as a function of orbital 
filling Tier. Right: Ground-state phase diagram for both HF and GW. The lines mark 
the instability for a transition from the paramagnetic (PM) to the ferromagnetic (FM) 
state. The orbital filling is n^ « 0.29 and t = \t';l]j (c.f., Ref. [55)1. 



For the variational ground state we can work with a wave function (^ that contains 
only diagonal parameters Ar,r- Non-diagonal parameters could only arise if we break 
the cubic symmetry or want to study states with magnetic orders not coUinear to the 
chosen spin-quantisation axis. Note, however, that for the study of spin excitations we 
must allow for non-diagonal variational parameters, see below. 

In our numerical analysis of this two-band model we will consider a tight-binding 
Hamiltonian Hq with generic hopping parameters which were already used in previous 
works and lead to the density of states at the Fermi energy shown in Fig. [1] (left). Due to 
the maximum in the density of states at approx. n^ = 0.29, in that range of band fillings 
there is the strongest tendency for a ferromagnetic state to be lower in energy than the 
paramagnet. This has already been demonstrated in Ref. [7]. Another important finding 
in that work is the huge importance of the exchange interaction J for the appearance of 
ferromagnetic order. This can be seen from the Gutzwiller phase diagram for our model 
in Fig. [1] (right). In contrast, the HF phase diagrams shows almost no dependence on 
the size of J; see also Ref. [H] where similar results have been reported for a two-band 
model in infinite dimensions. 



7.2. The magnetic susceptibility 

For the calculation of the spin susceptibility (I172p . we need to determine a Green's 
function matrix of the form (I178p in which the operators A^, i?^, B^ are given as 









hx 



B 



bi,b2,b{,b'2 



^k,{b2^fk+q,(ba) 



fl- 2-^^k ^k,{b>^^fk^q,{b\i)^ 
^' k 

_2 V^ b2,&i-t 

rr~ Z^ ^l^+<i ^k,{b'2\)^k+q,{b'^i) ' 



(193) 
(194) 
(195) 
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The matrix fll78p . which results from these operators is 4 + 16 + 16 = 36 dimensionaL 
Due to symmetries, this dimension can be reduced to 20 for a general wave-vector q. 
Along symmetry lines the symmetry reduction could even go further. In our numerical 
calculations, however, we did not exploit such symmetry considerations since the 
numerical efforts for a two-band model are still moderate, even in three dimensions. 

Note that there is a difference between Hartree-Fock and Gutzwiller RPA 
calculations concerning the elements of Il^{q,u) which have to be taken into account in 
our calculation of the susceptibility 

b,b' 

In Hartree-Fock RPA, due to the locality of the interaction terms in Hubbard models 
and the symmetries of Cg orbitals, the only elements of H°(g, w) which contribute are 
{{Alf^; {Al, fe/)^))[i, i.e., those which are diagonal with respect to the local orbital indices. 
This is different from the Gutzwiller RPA equations in which all Green's functions 
defined by the operators f ll93p - fll95p have to be taken into account. In particular. 
Green's functions ((^6^ ^2' ( 63 64)^))° with 61 7^ 62 or 63 ^ 64 cannot be discarded. 
The reason for this difference is the non-locality of the interaction matrix V^ in the 
time-dependent Gutzwiller theory. 

7.3. Results 

We prepare a ferromagnetic ground state in both HF and Gutzwiller approximation 
at band filling n„ ^ 0.2987 in order to be close to the maximum of the DOS. In 
general both schemes will give different magnetisations for the same set of interaction 
parameters. Therefore, one could either perform the comparison for fixed parameters or 
fixed magnetisation (cf. also Ref. [59]). To avoid this inconsistency, we present results 
for interaction parameters which lead to a fully polarised ferromagnetic ground state in 
both approximations, i.e. m = 2n„. Note that due to numerical reasons we have to stay 
slightly below this value in case of the Gutzwiller approximation. The corresponding 
interaction parameters are specified in the captions to Figs. 121 [3] which display the 
magnetic excitations obtained within both approximations. 

These spectra are composed of a low-energy magnon part due to the breaking 
of spin-rotational invariance and a high energy Stoner continuum which reflects the 
particle-hole spin-flip excitations of the 'bare' system, i.e H°(g, w), cf. Eq. f ll79p . For 
both methods we show the excitations along the (100) and (111) directions. The 
difference in these directions mainly arises due to the orientation dependence of the 
particle-hole dispersion which is significantly stronger along the diagonals. 

One first important difference between HF and Gutzwiller approximation concerns 

i~ I 

the difference in the magnetic band splitting A~ = Ep — Ep. In HF theory this value is 
just given by X~{HF) = {U -\- J)m and thus for strongly correlated systems produces a 
large gap 0{U) between the low energy magnon and a high energy Stoner continuum. 
On the other hand, we find \'{GA) is significantly reduced with regard to its HF 
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Figure 2. HF excitations for U = 10. Oi, J/U = 0.30 resulting in m = 0.5975 
(fully polarised) and A^ = Ep — Ep w 7.771. The magnon dispersion is fitted by 
emaanon{q)/t = D\(f\'' {I + p\(f\'') witli D^""^ « D^^^ = 100 X 10-3^ Scaling: t = \t[ 
(c.f., Ref. [58]) and qx,y,z G (— ti','?''). 



''ddcr\ 



counterpart. For the present system \^{GA) k. 1/AX~{HF). Given the broadening of 
the Stoner continuum with increasing transferred momentum, the low energy magnon 
thus rapidly merges with the continuum in the time- dependent GA as can be seen from 
Fig. [HI As a consequence the excitation at a; ~ J, corresponding to a respective spin-flip 
in the two orbitals, is only visible in HF-I-RPA along the (100) direction, whereas in 
the time-dependent Gutzwiller approach it is already within the continuum. Note that 
the overestimation of the Stoner excitation energy within HF+RPA is a longstanding 
problem in solid state theory as discussed in Ref. [60] . 

At g* = all the weight is contained in the zero frequency Goldstone mode. 
The existence of this excitation provides an important consistency check of the 
Gutzwiller+RPA approach similar to the analogous finding in HF+RPA. The positive 
dispersion of the magnon further demonstrates that the underlying Gutzwiller solution 
is a stable energy minimum which is not destroyed by the fluctuations. The spin-wave 
stiffness, i.e. the quadratric coefficient of the magnon dispersion, is significantly larger 
in HF+RPA than in time-dependent Gutzwiller theory. Note, however, that to a certain 
extend this huge difference is caused by an instability of the ferromagnetic ground state 
with respect to an incommensurate phase which is found for interaction parameters not 
much smaller than those used in Fig. [3] 
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\q = (qqq) 



a/t 



Ql/t 



Figure 3. GW excitations for U = 10. Ot, J/U = 0.30 resulting in m = 0.5728 (almost 
fully polarized) and A^ = Ep — Ep sa 2.128i. The Magnon dispersion is fitted by 
emagnon{q)/t = D\^'^ {I + I3\q\'^) with D^^^ = 1.34 X lO'^ and I?!" ^ 1.30 x 10"^ for 
small wave vectors. Scahng: t = [t^j^^l (c.f., Ref. [58]) and qx,y,z G (— 7r,7r). 



8. Summary 

In this paper we have given a detailed derivation of the time- dependent Gutzwiller 
approximation for multi-band Hubbard models. The basic assumptions which underlie 
the method can be summarised as follows. First, it is assumed that the dynamics 
of the Slater-determinant (upon which the Gutzwiller projector acts in the starting 
Ansatz) is determined by the so-called Gutzwiller Hamiltonian, Eq. (I140p . which leads 
to an equation of motion similar to standard RPA, Eq. f ll38p . Second, the dynamics 
of the variational parameters is determined from the assumption that at each instant 
of time the energy is minimised. This leads to a linear relation between variational 
parameter and density fluctuations, Eq. (I154p . Third, as in the standard HF+RPA 
approach it is assumed that the external perturbation and thus the density fluctuations 
are small, Eqs. fll06p . fll07p . We have seen that also in the multi-band case these 
assumptions lead to a consistent theory in the sense, that an instability which is signalled 
within the Gutzwiller+RPA corresponds to a (second-order) phase transition which one 
would obtain from the bare variational Gutzwiller approximation. We have further 
demonstrated that for ferromagnetic ground states the Gutzwiller-fRFA leads to the 
appearance of the Goldstone mode as expected for systems which break continuous spin 
symmetry. 

The formalism as developed in its present form can now be straightforwardly applied 
to the investigation of correlation functions in strongly correlated multi-band systems 
as e.g. pnictides, manganites, cobaltates etc. On the other hand, a natural application 
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of the theory would also comprise the investigation of e.g. orbital quenches for which 
the small amplitude assumption for the density matrices has to be abandoned. For 
single-band Hubbard models such a fully time-dependent formulation of the Gutzwiller 
approximation has been recently presented by Schiro and Fabrizio [611 [62] where also 
the second assumption above has been replaced by separate equations of motion for 
the variational parameters. Future work should thus address the question whether 
their approach reduces to the present theory in the small amplitude limit and how it 
eventually can be extended to the multi-band case. 

Appendix A. Second order expansion of determinants 

According to Eq. (iT6l) . section 13.11 expectation values m'j p can be written as 
determinants of certain matrices A with elements that are linear functions of the local 
density matrix C^^y. In the variational ground state, C°,^/ is diagonal, and, if we chose 
a proper order of the orbitals 7, to 0-th order the matrix A is diagonal too, 

Al^ ... 

0' Al, ... 



A%j^ 



m^r r, = \A\ = \A° 



"ij 



(A.l) 



For J 7^ /' at least one of the diagonal elements A^^ vanishes and we find 

i 76/ yG{l...N)\I 

as expected. In order to calculate the first and second derivative of m^ p we need to 
expand the determinant 

\A\ = \A^ + 6a\ (A.3) 

up to second order with respect to the matrix elements Saij. For this expansion one 
readily finds 

\A\ — \A^\ = |y4°l y^ ''' _L mOi y^ Saj^jdajj -\- 6aij6aj^i 



«j 



«,« jj 



Note that for A'^^ = the right-hand side is defined by the corresponding limit A^- — ^ 0. 

Appendix B. Invariance of Second-Order Expansions 

Appendix B.l. Equivalence of the Lagrange- functional expansion 

In this section, we show that the interaction kernel Kyyi i^i (H64p is identical to Kyyi 
in Eqs. f ll55p - fll56p . To this end, we choose again some arbitrary independent and 
dependent variational parameters A^ and A^, c.f. Eq. f ll23p . By construction, the 
constraints (I122p are automatically fulfilled as a function of A^ and p, i.e., we have 

gn{\\\\p),\\p) = Q. (B.l) 
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Consequently, all first or higher-order derivatives of (IB.ip with respect to A^ and py 
vanish. For example, the first-order derivatives lead to 

dgn dgn 



d\\ 
dgr. 



d\\ 






dpY dpY 
Using the matrices 

_ dgn 



dgn ^ dgn d\% 



G 



n,X 



9Ai' 



X 



R 



d\% dpY 



dX% 



x,z 



X 



dX\ 



Qx, 



dXi 



Y 



X 



dpY 



we can write Eqs. (IB.2p -f lBl3|) as 
dgn 



w,-' f^^]-^ ' 



dgr, 



- [GQl 



(B.2) 
(B.3) 

(B.4) 

(B.5) 
(B.6) 



dpY L-~.Jn,y 

With the classification of dependent and independent variables we are in the 
position to evaluate the antiadiabaticity conditions fll62p - fll63p . First, Eq. f ll62p leads 
to 

dgn 



Eit«uEir*Ak+E|r*'' 



X S^-V 



Y 







Y 



which, together with Eqs. (lB.5p - (lB.6p . yields 



a 



SX'^ - R6X' - Q6p 



0. 



(B.7) 



(B.8) 



Since the square matrix G should be invertible, the bracket in (IB.SP must vanish. Hence, 
we find the relation 

^A'^ = R6X + Q6p (B.9) 

which determines the dependent-parameters fluctuations ^A'^ as a function of (5A' and 
6p. 

Applying the separation of dependent and independent parameter fluctuations to 
the second set of Eqs. f ll63p yields 

G^ Ld' L^d 



<5A' 
5A^ 



UP 



^Ap 



5p 



(B.IO) 



Here we introduced the six matrices 



r" — 
^z,z' — 



' dp 



, . . . , i^xy 



d^L 



(B.ll) 



dX'zdX'z, ""'^ dX%dpY 

of second derivatives. With ( 1B.9P and the second 'row' of Eqs. (IB.lOp one can write the 
Lagrange-parameter fluctuations as a function of 6X^ and 6p, 



5 A = -[G^]-^ \{L'^' + L'^'^R)6X' + {L'^'^Q + L'^P)6p 



(B.12) 
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we eventually find 
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(B.13) 



6X' 






(B.14) 



Equations (]B.14p . f lB.12p . and (JB.9P now enable us to write all fluctuations 6\\ 6X'^, 
and 6A as functions of the density fluctuations 6p. These relations can be inserted into 
the second-order expansion of the Lagrange functional, 



+ 



{{5pfLP''5X'' + {SpfU'SX + {5XYL"'5X') 



+ 



+ 2{5KYG 



5X'^ - R5X' - Q5p 



(B.15) 



in order to calculate Kyy' iii Eq. fll64p . However, to prove just the identity of Kyy' 
and Kyyi iii (H55p it is sufficient to apply only Eq. (]B.9p to the expansion ( jB.lSp . This 
leads to 



2(5L(2) = {6pf{LPP + Q'^L'^P + LP'^Q + Q^L''''Q)6p 

+ {6XY{D' + U'^R + FfL''' + lfL''^R)6X 



(B.16) 



+ 



\T 



{5pf{LP' + LP'R + g^L^' + Q^L'^'^R)6X' + 

As we will show below, the matrices fll48p -f fT50|) which determine the second order 
expansion (115ip are the same as the corresponding matrices in ( IB.lGp . Hence, we have 

5^(2) = 5L(2) . (B.17) 

Since the antiadiabaticity condition 



d5X\ 



dSX\ 







(B.18) 



for 5E^'^^ reproduces Eq. flB.14p . the identity of Kyy ^-^id Kyy is then finally 
demonstrated. 

It remains to be shown that the matrices fll48p -( 1T50|) agree with those in (]B.16p . To 
this end, we use the explicit form f ll24p of the energy functional f ll20p that appears in 
the definition of the matrices fll48p - fll50p . As an example, we consider the matrix 
iVf'' and show that it is identical to the matrix in the first line of flB.16p . With 
similar derivations one can prove the same for the other matrices f ll49p .( TT50l) and their 
counterparts in (lB.16p . 

Using flT24l) and f lT48l) we find 



M^'y, = [E^P + Q^E^P + W'^Q + Q'^E'"'Q]yy, + 2 V |-p 



d^Xj, 



dpydpY' 

(B.19) 



The time- dependent Gutzwiller theory for multi-hand Hubbard models 
Here, the matrices 



~dd 



32 



(B.20) 



and g^'^, . . . , g^'^ are defined as in (IB. Ill) only with L replaced by E or gn respectively. 
Obviously, the matrix in the first line of flB.lOP is identical to M^^ if 



X 



dXx dpydpY' 



To prove (IB.2ip . we use the fact that the second (total) derivatives of (IB.1|) with respect 
to the densities py vanish 

dgn 



dpydp 



Y> 



V€ + (r~9i' + ~<'Q + Q^~9'.'Q]y,y, + 2 E 1?^ ^'^^ 



X 



dXx dpydpyi 







Equation ( IB. 211) is therefore fulfilled if 



X ^^^i 



d'X'x 



0. 



dXxJ dpydpyi 

This equation, however, holds trivially, since fll58p leads to 

dgr, 



dL dE 



+ Ea- 







d\z dXz ^-^ d\z 

n 

for all parameters \z and in particular for \z = Xx as it appears in (]B.24|) . 



(B.22) 



(B.23) 



(B.24) 



Appendix B.2. Linear transformations of the density matrix 

In investigations of our translationally invariant lattice systems ([1]) it turns out to be 
more convenient to work with fluctuations 5p which are linearly related to the density- 
matrix fluctuations, 

6p = E-6jl (B.25) 

c.f., Eqs. (lC.27p - (IC.28p and the resulting Green's functions f ll78p . The effective second- 
order functional fll55p - fll56l) in terms of the fluctuations 6fl is then given as 



6E<^'\6fl) = l{6flf{E^k^^E)6p 



(B.26) 



with K'^'^ as deflned in (11561) . For numerical calculations it is important to show that 
one obtains the same kernel 



K^^ = 'E'^ K^'' 



(B.27) 



as in (JB.26P if the transformation (JB.25P and the antiadiabaticity condition are applied 
in the reverse order: If we apply ( JB.25P flrst to (I15ip . we obtain 



2 



{6flfE^M^^E5fI + 2{6XYM^f'E6fI+{6XYM^^6X' . (B.28 
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The antiadiabaticity condition for Sjl then reads 



6X' 



M 



XX 



M^PESfl . (B.29) 



Inserted into ( JB.28P this equation yields 



6E^^\6fI) = Eo + ^{6flfK^^SfI , (B.30) 



M 



XX 



M^^H = E'^RPpE 



as claimed above. 



Appendix C. Explicit form of the second-order expansion 

We calculate the second-order expansion of the Lagrange functional with respect to the 
variational parameters Aj;r,r' and the density matrix (!65|) . For the general consideration 



in section [5] and Appendix B it was convenient to subsume the parameters Ar,r' and their 
conjugates App, in a set of n^ parameters \z, c.f., Eq. f ll2ip . Here in this appendix, 
where we aim to resolve the explicit structure of the second-order expansion, it is better 
to take the difference between Ar,r' and App, into account. 

The constraints ( 135|) -( 136|) . the local energy (12^ . and the renormalisation matrix (!27|) 
are all functions only of A*.p p/ , Aj;r,r' and of the local density matrix C°^ ^, . For simplicity 
we use the joint variables A^, {AD* for all these local variables, i.e., it is either 

Al = Al^^^^ = {cl^c.^J or a; = A^,p, = A,r,r' . (C.l) 

With respect to the parameters A*.pp,, Aj;r,r' the second derivatives of (I35|) - fl36|) . f l24|) . 
and ( 1271) are quadratic functions of the form ~ {Al)*Al,. Due to the Hermiticity of the 
density matrix the same can be achieved with respect to the local density matrix. Then 
the only finite second derivatives of the Lagrange functional 

L = T + J2 E^,U{iAly}, {Al}) + J2 ^^,ng^A{iAlr}, {A^}) (C.2) 

i i,n 

= T -\- Lioc 



are 



i¥=j 



^ ^0 (C.4) 



whereas 






d{Aiyd{Ai,)* dAldAi, °" ^^-^^ 

The second-order expansion of the constraints and the local energy is 
straightforward since only local fluctuations SA^ couple, 

^A2 = EE(^^')*^S"^'^^'' (C.6) 
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T^IOC 



d^L 



loc 



and the Fourier transforms of the local fluctuations 
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(C.7) 



{C.i 



All derivatives in this section (e.g., ( 1C.7P ) have to be evaluated for the ground-state 
values of the variational parameters Aj;r,r'; the density matrix p, and the Lagrange 
parameters Aj^„. Note that the density-matrix fluctuations SA^^ ^^ can be written as 



^K... 



^E^"^-^(C.^...) = ^E^(4 



where the operator A^ has been defined in (11691) . 



.c 



k,ai k+q,a2 



(C.9) 



In addition to flC.6p . we need to take into account the mixed terms ~ 5y4^(5Aj„. In 
real space, their contribution is given as 



«f = E ( 






(CIO) 



If we introduce the Fourier transforms 5A^ of the fluctuations 5Aj_„, we can write (IC.lOll 
as 

q n,v 

Here, we used that the constraints gi^n are assumed to be real and lattice-site independent 
such that 

More involved than the calculation of (1C.6I) is the expansion of the kinetic energy. 
Here we find 

with 



(C,13) 



*7^i 0-1,0-2, 0-^,0-2 ^y 

1 / 5<ii ^ (C 






S(AVrSAi, ^"'i^l ^'^'>'''^' 



+ - 



2 \diAiy dAl, 

+ c.c. 



<^2i dq'^i ^Uj>2, 



Mt d{AlY 



(C.14) 
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i¥=j 0-1,0-2,0-^,0-^ 
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(C.15) 



E 



''-^ ^<U*(^^^)*+<i,^M^(^^i) 



+ c.c. . 



The fact that the complex conjugates give the terms not exphcitly shown in Eqs. flC.141) - 
flC.lSp follows from the relations 



d{A,ydA,, 



d{A,y ' 
diA,>)*dA, 









(C.16) 

(C.17) 
(C.18) 



For our translationally invariant ground state it is more convenient to write (10.14^ - 
( JC.151) in momentum space: With the Fourier transforms of the local fluctuations the 
term (]C.14p reads 

^^/'^ = E E(^^')K-,-' + «.',J*]K' (c-19) 

q v,v' 



where 



K. 



1 

q;v,v' 



E 






d\l\ 



^^"^'"^'"'^'"^5(A)*M,. 



c 



(C.20) 



Here we assumed that the renormalisation matrix is lattice-site independent and 
introduced the tensor 



with 



^oi,02,ol,o^(g) = 7-E^fc+r(4oi^feo^) 



0-1,02 _ \ ^^oi,02 ifc(-R;--R,) 



(C.21) 



(C.22) 



Note that for g = the tensor (!C:2T1) . 

-Eoi,o2,o;,o^ = -^^1,02,0^,0^(0) (C.23) 

has already been defined in ( !32l) . For the evaluation of the second ('transitive') term 
(1C.15I1 we write the non-local density-matrix fluctuations as 



^(4iS.,) = fE^^^''-'-'^"^^(4oi%.,). 



(C.24) 



k,k' 
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(C.25) 



q,fc v;a[,a'2 



with 






r dqa\ 



^<T2 ) ^k+q 



d (qll 



<y\ ,0-2 I & 



.o'l,o"2 



d{A,y 



(C.26) 



In principle, Eqs. f lC.25|) -( IC.26|) allow us to calculate all second-order couplings 
of density-matrix and parameter fluctuations that arise from 5T^ . For numerical 
calculations, however, these equations are not very useful due to the explicit 
k dependence of flC.26p . It is much easier to introduce the two auxiliary fluctuations 



631 = 63" , , 

^ cr2,cri,(T2,cr-j^ 



631=63" , , 

W (72, Ui, 0-2,01 



^' k 
T^Z-^'^k+q " \'^ka['^k+qa'J 



(C.27) 
(C.28) 



where w = {a2,<yi,o'2,o'[) is an abbreviation for quadruples of indices a. With these 
deflnitions we can write (IC.25P as 



^^f ^ = EE \i5AirKl[^^63l + {6AiyKl^^)6B: 



(C.29) 



q v,w 



where 



K 



t(i) 



+ {63iy{Kl^^^r6Al + {6Biy{Kl^^^y6Al 
, d (qol 



= €]■ 



v(o2,oi,a'2,o[) "ioi Q(j\\* 



K 



t(2) 



dql\ 



(C.30) 



(C.31) 



Note that we introduced the two different fluctuations ( ]C.27|) . (1C.28P only because they 
allow us to write the second-order expansion in a relatively simple form. In fact, these 
fluctuations are not independent but related through 

63" , , =(63-" , ,y . (C.32) 

01,02,Oj^,02 ^ 02,0\,02,0^' \ ' 

Altogether we end up with the following second-order expansion of the Lagrange 
functional 

/ 6A'i \ 



6L^'^ = j-Y,{6A'^ 63^ 6W 6A" )* K" 



63'i 

63" 

\6l" J 



(C.33) 
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where 



K'^ = 
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(C.34) 



(C.35) 
(C.36) 
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As described in section I5.4[ the antiadiabaticity condition leads to an effective second- 
order functional only of the density matrix. This condition can be evaluated directly 
for the second-order expansion flC.331) since the fluctuations SA"^, 6B'^, 6B are some 
linear functions of the density-matrix fluctuations 5{cl,ai'^k+q,a2)^ '^•^•' Appendix B.2 To 
this end, we distinguish the fluctuations of the local density matrix SA'j^ and of the 
variational parameters SA^ as well as the corresponding blocks in the matrix flC.34p . 
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The resulting functional is then given as 
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where 
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Note that V^ (unlike K'^) includes finite couplings also between the fluctuations 6B'^, 
6B . The calculation of V^ (for fixed q) only involves the handling of finite-dimensional 
matrices. In contrast, the evaluation of the functional (10.25^ (instead of (JC.291) ) would 
have lead to significantly more complicated equations. 
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Appendix D. Explicit form of the Gutzwiller-RPA equations 

In this appendix, we prove that the general Gutzwiller RPA Eqs. fll44p lead to the 
Green's function matrix fll79p if applied to our multi-band Hamiltonian ([1]). With the 
abbreviations SD"^, D1^ for the three fluctuations 6AI, SB'^, SB^ and the corresponding 
operators A^, 5^, B^, we have to show that the Green's function matrix 

U,,,,{q,u) = {{bl;{bl,y))^, (D.l) 

as given in (I179p . obeys the equation 

sDi = J2((^i;{bi.y))jfy (D.2) 

Using the explicit form fll79p of Il{q,u), this equation can also be written as 

Note that the excitation amplitudes Sf^ enter the problem through the perturbation 
operator 



5Vf 

cri, (72, cr'^,cr'^ 



k ^, ^„ ^' ^' 



which is needed to deflne the general Green's functions (11780 . 

Before we prove Eq. (ID.Sp . it is instructive to consider the case \^'' = in which 
the three fluctuations 6AI, SB'^, 6B!l, are decoupled and we can set f^''^ = /^'"^ = 0. 
We start this derivation in the eigenbasis of the Gutzwiller Hamiltonian fll80p . It leads 
to the simplest form of the matrix E in Eq. f ll44p which then reads 

(c - (E,+,,„, - i?fc,.j)5(/.^^ V,,,^)^P/P^ (D.5) 

-V''k,a2 ^ ^k+q,ai)^t{k+q,a-i),{k,a2) ■ 



Here the excitation amplitude is given as 

Sf^k+q,a,Uk,a2) = E ^fZMuJ*<a2 ■ (D-6) 



0'l,o-2 



Note that the factor n^^^^ — n^k+q,ax — '^^ i'^ ( ID.SP represents the particle-hole and the 
hole-particle channels in Eq. (I144p . For simplicity, we will drop the corresponding labels 
hp/ph in the following. 

With the transformations flT82D . flT83D . Eq. flD3|) leads to 
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k '*i'"2 



As expected, we therefore find 



^^^..2 = E((^^.^.'Ki,.^)'))°^A,4 ^D-8) 



with the ('retarded') Green's function 
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as introduced in (11841) . 

Now we consider the case of a finite interaction matrix V. Using our abbreviation 
bDy, for the amphtudes bA^, SB^,, 53^ the Lagrange functional 5L^^) has the form 

6L^'^=Y.(^Dirv;!yiSDl,). (D.IO) 

q,fi,fi' 

With this additional interaction term, Eq. fID.SP reads 



(c. - (E,+,,„, - i5.,„j)5(/il,,^/i,^^,„^) + {nl^^ - nl^^^J (D.ll) 
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where 
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and 



"^ "■'(Ti,(T2,CTi,4''fe "T "iCTi,a2,<,4 fc+g 
The derivatives in (1D.13P can be further evaluated using the transformation (I182l) . (ll83p . 
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Depending on the particular fluctuations SDj^, the remaining derivatives on the r. h. s. 
of equations (ID.15|) . flD.16|) are given as 
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With Eqs. ( 1D.12p -( rD.19P we are now in the position to evaluate ( ]D.11|) . To this end, we 
proceed as in (JD.71) . 
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The sums over /i and /i' lead to nine contributions which can all be evaluated using Eqs. 
flDT7l) -f lD19|) . As a result we find 

m + Y.^{Al- {D^^'))lVl^,6Dl, = J2{{Al; {DlV)Z SP, . (D.21) 
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where the 'non- interacting' Green's functions {{Al; (-D?^)^))[i in flD.2ip are given as in 
(ID.gp . apart from additional factors e^^'*^" or el 

U^^oi.oa' (-°03,04, 0^,04) //ui 

{{Al^^^^;{Bl^^,,y)Z 
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With flD.2ip . we have proven the 'first' set of Eqs. fID.Sp . i.e., those with fj, = v 
If we replace 6A„^^„^ in the first line of Eq. ( ]D.20p by 
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or by 
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the remaining Eqs. (1D.3P are derived in the very same way as (lD.2ip . This closes our 

proof of Eq. flT79|l . 
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